0. STUDY CONTEXT AND OBJECTIVES

Tile: Telepressure & Suicide in Secundary School Spain - Network Analysis: A Network Analysis Approach

Authors:
José-Manuel Guerra¹,
Jesús Martínez-Calvo²,
Diego Diaz-Milanes^2,3

Affiliation:

¹ Department of Social Psychology, University of Seville, Seville, 41018, Spain

² Department of Quantitative Methods, Universidad Loyola Andalucía, Sevilla, Spain

³ Health Research Institute, University of Canberra, Canberra, ACT, Australia

Reproducibility Statement

This R Markdown file contains the de-identified data and all R code used in the study, ensuring full transparency and reproducibility.

Analysis Workflow

The workflow comprises ten steps:

Data import and preprocessing
Load the dataset, screen missingness/outliers, and prepare variables.
Descriptive statistics
Summarize sociodemographics, inspect the distribution of variables, and explore gender differences where relevant.
Undirected network modeling
Estimate a partial correlation (GGM) network of WTM items using EBICglasso.
Network invariance testing
Compare network structure, global strength, and specific edges across gender groups via the Network Comparison Test (FDR-adjusted edgewise tests).
Directed network modeling
Learn a Bayesian Network (DAG) for potential directional dependencies among WTM items with bootstrap aggregation.
Predictive evaluation
Assess model predictability (e.g., \(R^2\), RMSE) via cross-validation and compare performance across models.

1. DATA LOADING & PROCESSING

In this section, the packages, function design, and data used are listed. Then, data processing—including variable selection and filtering by units of analysis—was performed.

1.1. Import Packages & Functions

# ---- Load libraries (already attached above; keep here if you knit this chunk standalone)
suppressPackageStartupMessages({
  library(haven); library(tidyr); library(dplyr); library(corrplot); library(psych)
  library(parameters); library(ggplot2); library(semTools)
  library(tibble); library(purrr); library(effectsize); library(EGAnet); library(qgraph)
  library(NetworkComparisonTest); library(bootnet); library(mgm); library(bnlearn); library(caret)
  library(forcats); library(knitr); library(kableExtra)
})

# ---- Helper: safe version + loaded flag
safe_version <- function(pkg) {
  v <- tryCatch(as.character(utils::packageVersion(pkg)), error = function(e) NA_character_)
  if (is.na(v)) "" else v
}
is_loaded <- function(pkg) pkg %in% loadedNamespaces()

# ---- Packages + key functions actually used/reported
pkgs <- tibble::tibble(
  Package   = c(
    "readxl","tidyr","dplyr","psych","semTools","ggplot2","effectsize",
    "EGAnet","qgraph","NetworkComparisonTest","bootnet","mgm","bnlearn",
    "caret","forcats","knitr","kableExtra","lavaan","semPlot","corrplot","haven","parameters"
  ),
  Functions = c(
    "`read_excel`",
    "`pivot_longer`, `pivot_wider`, reshaping support",
    "`select`, `mutate`, `filter`, `group_by`, `summarise`, `bind_rows`",
    "`describe`, `polychoric`, `alpha`, `omega`, `splitHalf`",
    "`mardiaSkew`, `mardiaKurtosis`, invariance helpers",
    "`ggplot`, `aes`, `geom_line`, `labs`, `theme_bw`",
    "`cohens_d`",
    "`bootEGA`, `dimensionStability`",
    "`qgraph`, `EBICglasso`",
    "`NCT`, `plot`",
    "`estimateNetwork`, `bootnet`, `corStability`",
    "`mgm`, `predict`",
    "`boot.strength`, `averaged.network`, `bn.fit`, `predict`",
    "`confusionMatrix`",
    "`fct_reorder`",
    "`kable`",
    "`kable_styling`, `add_header_above`",
    "`cfa`, `fitMeasures`",
    "`semPaths`",
    "`corrplot`",
    "`read_sav`, `as_factor`",
    "`model_parameters`"
  )
)

# ---- Add Version and Loaded columns
packages_info <- pkgs %>%
  mutate(
    Version = vapply(Package, safe_version, character(1)),
    Loaded  = ifelse(vapply(Package, is_loaded, logical(1)), "Yes", "No")
  ) %>%
  select(Package, Version, Loaded, Functions)

# ---- Display table
knitr::kable(
  packages_info,
  caption = "Table 1. Packages used, versions, loaded status, and key functions"
) %>%
  kableExtra::kable_styling(full_width = FALSE, position = "center")

Table 1. Packages used, versions, loaded status, and key functions
Package	Version	Loaded	Functions
readxl	1.4.3	No	`read_excel`
tidyr	1.3.1	Yes	`pivot_longer`, `pivot_wider`, reshaping support
dplyr	1.1.4	Yes	`select`, `mutate`, `filter`, `group_by`, `summarise`, `bind_rows`
psych	2.4.1	Yes	`describe`, `polychoric`, `alpha`, `omega`, `splitHalf`
semTools	0.5.6	Yes	`mardiaSkew`, `mardiaKurtosis`, invariance helpers
ggplot2	3.5.2	Yes	`ggplot`, `aes`, `geom_line`, `labs`, `theme_bw`
effectsize	0.8.9	Yes	`cohens_d`
EGAnet	2.0.4	Yes	`bootEGA`, `dimensionStability`
qgraph	1.9.8	Yes	`qgraph`, `EBICglasso`
NetworkComparisonTest	2.2.2	Yes	`NCT`, `plot`
bootnet	1.5.6	Yes	`estimateNetwork`, `bootnet`, `corStability`
mgm	1.2.14	Yes	`mgm`, `predict`
bnlearn	4.9	Yes	`boot.strength`, `averaged.network`, `bn.fit`, `predict`
caret	6.0.94	Yes	`confusionMatrix`
forcats	1.0.0	Yes	`fct_reorder`
knitr	1.45	Yes	`kable`
kableExtra	1.4.0	Yes	`kable_styling`, `add_header_above`
lavaan	0.6.17	Yes	`cfa`, `fitMeasures`
semPlot	1.1.6	No	`semPaths`
corrplot	0.92	Yes	`corrplot`
haven	2.5.4	Yes	`read_sav`, `as_factor`
parameters	0.22.1	Yes	`model_parameters`

Estamos utilizando 22 paquetes.

Design and generation of a formula for plotting centrality measures of two networks.

compareCentrality <- function(net1, net2,
                              include = c("Strength",
                                          "Closeness",
                                          "Betweenness",
                                          "ExpectedInfluence",
                                          "all",
                                          "All"),
                              orderBy = c("Strength",
                                          "Closeness",
                                          "Betweenness",
                                          "ExpectedInfluence"),
                              decreasing = T,
                              legendName = '',
                              net1Name = 'Network 1',
                              net2Name = 'Network 2'){
  
  if(include == "All" | include == "all"){
    include = c("Strength",
                "Closeness",
                "Betweenness",
                "ExpectedInfluence")
  }
  
  df <- centralityTable(net1, net2) %>% filter(measure %in% include)
  
  df %>% 
    mutate(graph = case_when(graph == 'graph 1' ~ net1Name,
                             graph == 'graph 2' ~ net2Name),
           graph = as.factor(graph),
           node = as.factor(node)) %>% 
    
    mutate(node = fct_reorder(node, value)) %>% 
    
    ggplot(aes(x = node, y = value, group = graph)) +
    
    geom_line(aes(linetype = graph), size = 1) +
    
    labs(x = '', y = '') +
    
    scale_linetype_discrete(name = legendName) +
    
    coord_flip() +
    
    facet_grid(~measure) +
    
    theme_bw()
  
}

1.2. Data Loading

The dataset from the #### project (Ethical Committee: ####) was loaded, including only the minimum set of variables required for the present study.

df_general <- read_sav("Estudio Telepresión - Fase 3 - ESO.sav")
df_general

1.3. Filtering

Filters were applied to remove participants who did not complete the isntrument or indicate their gender, and to exclude cases with missing values:

tele_df <- dplyr::select(df_general,c(Sexo,EdadV1,V2,AUTO:TPB))
names(tele_df)[1:3] <- c("sex","age","academic_year")
tele_df <- as.data.frame(tele_df)

tele_df$sex <- as.factor(case_when(tele_df$sex == 1 ~ "Male",
                         tele_df$sex == 2 ~ "Female",
                         tele_df$sex == 3 ~ "I do not know"))

tele_df$academic_year <- as.factor(case_when(tele_df$academic_year == 1 ~ "First academic period",
                         tele_df$academic_year == 2 ~ "First academic period",
                         tele_df$academic_year == 3 ~ "Second academic period",
                         tele_df$academic_year == 4 ~ "Second academic period"))

tele_df$PS <- as.factor(case_when(tele_df$PS == 0 ~ "Passive-Agressive",
                         tele_df$PS == 1 ~ "Assertive",
                         tele_df$PS == 2 ~ "Agressive",
                         tele_df$PS == 3 ~ "Passive"))

tele_df$TPB <- as.factor(case_when(tele_df$TPB == 1 ~ "Low risk",
                         tele_df$TPB == 2 ~ "Medium risk",
                         tele_df$TPB == 3 ~ "High risk"))

tele_df <- dplyr::select(tele_df, c("sex","academic_year","PS","TPB","age","AUTO","HETERO",
                                    "TP","SP","RS","CIUS","IDS","AB","VB","VCB","ACB","AF",
                                    "AA","APS"))

initial_rows <- nrow(tele_df)
tele_df <- na.omit(tele_df)

print(paste("Initial number of rows:", initial_rows, "| Final number of rows:", nrow(tele_df)))

## [1] "Initial number of rows: 820 | Final number of rows: 813"

tele_df

1.4. Variable Selection

Creation of dataframes tailored to meet the requirements of the forthcoming data analyses.

1. Dataframe for the analysis of individual items

The following table displays the first five rows of the set of items from the WTM instrument, used in the present study.

# Prepare the dataframe
network_variables <- dplyr::select(tele_df, age:APS)
pieColor <- rep("#EB9446", length(network_variables))

# Display first 5 rows with a formatted table
head(network_variables, 5) %>%
  kable(caption = "Table 2. Set of items from the WTM") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 2. Set of items from the WTM
age	AUTO	HETERO	TP	SP	RS	CIUS	IDS	AB	VB	VCB	ACB	AF	AA	APS
14	10	27	3.500000	2.75	3.333333	2.533333	2.4	1.333333	1.333333	1.000000	1.000000	2.25	3.00	2.75
13	55	7	2.500000	2.25	2.000000	2.866667	3.2	1.666667	2.333333	1.000000	2.000000	1.00	1.25	1.25
13	55	53	1.833333	3.25	1.666667	1.333333	1.2	1.000000	2.000000	1.666667	1.000000	4.00	1.50	4.00
14	50	53	2.333333	2.75	2.000000	1.800000	1.0	1.666667	2.000000	1.333333	1.333333	4.00	4.00	3.75
14	45	47	2.500000	2.00	1.666667	2.600000	3.8	1.666667	1.333333	1.000000	1.000000	1.25	2.75	2.75

2. Dataframe for exploring sex differences

The following table shows the first five rows of the dataset including the WTM item responses, the total score, and participants’ sex

# Select CEI items, total score, and gender
NV_sex <- dplyr::select(tele_df, c(sex,age:APS))

# Display first 5 rows with styling
head(NV_sex, 5) %>%
  kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 3. WTM items, total score, and gender variable
sex	age	AUTO	HETERO	TP	SP	RS	CIUS	IDS	AB	VB	VCB	ACB	AF	AA	APS
Female	14	10	27	3.500000	2.75	3.333333	2.533333	2.4	1.333333	1.333333	1.000000	1.000000	2.25	3.00	2.75
Female	13	55	7	2.500000	2.25	2.000000	2.866667	3.2	1.666667	2.333333	1.000000	2.000000	1.00	1.25	1.25
Female	13	55	53	1.833333	3.25	1.666667	1.333333	1.2	1.000000	2.000000	1.666667	1.000000	4.00	1.50	4.00
Male	14	50	53	2.333333	2.75	2.000000	1.800000	1.0	1.666667	2.000000	1.333333	1.333333	4.00	4.00	3.75
Female	14	45	47	2.500000	2.00	1.666667	2.600000	3.8	1.666667	1.333333	1.000000	1.000000	1.25	2.75	2.75

2. Dataframe for exploring academic year differences

The following table shows the first five rows of the dataset including the WTM item responses, the total score, and participants’ academic year

NV_academic_year <- dplyr::select(tele_df, c(academic_year,age:APS))

# Display first 5 rows with styling
head(NV_academic_year, 5) %>%
  kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 3. WTM items, total score, and gender variable
academic_year	age	AUTO	HETERO	TP	SP	RS	CIUS	IDS	AB	VB	VCB	ACB	AF	AA	APS
First academic period	14	10	27	3.500000	2.75	3.333333	2.533333	2.4	1.333333	1.333333	1.000000	1.000000	2.25	3.00	2.75
First academic period	13	55	7	2.500000	2.25	2.000000	2.866667	3.2	1.666667	2.333333	1.000000	2.000000	1.00	1.25	1.25
First academic period	13	55	53	1.833333	3.25	1.666667	1.333333	1.2	1.000000	2.000000	1.666667	1.000000	4.00	1.50	4.00
First academic period	14	50	53	2.333333	2.75	2.000000	1.800000	1.0	1.666667	2.000000	1.333333	1.333333	4.00	4.00	3.75
First academic period	14	45	47	2.500000	2.00	1.666667	2.600000	3.8	1.666667	1.333333	1.000000	1.000000	1.25	2.75	2.75

2. DESCRIPTIVE ANALYSIS

2.1. Sociodemographic

Participants distribution by sex and age:

# sex distribution
sex_table <- table(tele_df$sex)
sex_percent <- round(prop.table(sex_table) * 100, 2)

sex_df <- data.frame(
  sex = names(sex_table),
  Frequency = round(as.numeric(sex_table), 2),
  Percentage = paste0(sex_percent, "%")
)

# Age summary
df_age <- tele_df[!is.na(tele_df$age), ]
age_summary <- summary(df_age$age)
age_sd <- round(sd(df_age$age), 2)

age_df <- data.frame(
  Statistic = names(age_summary),
  Value = round(as.numeric(age_summary), 2)
) %>%
  bind_rows(data.frame(Statistic = "Standard Deviation", Value = age_sd))

# Display tables
kable(sex_df, caption = "Table 5. sex distribution (frequency and %)") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 5. sex distribution (frequency and %)
sex	Frequency	Percentage
Female	400	49.2%
I do not know	10	1.23%
Male	403	49.57%

kable(age_df, caption = "Table 6. Age summary statistics") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 6. Age summary statistics
Statistic	Value
Min.	12.00
1st Qu.	13.00
Median	14.00
Mean	13.84
3rd Qu.	15.00
Max.	18.00
Standard Deviation	1.24

2.2. Variables

Univariate descriptive analysis of the virables included in the analysis:

# Descriptive statistics
descrgen <- describe(network_variables)

# Build dataframe with selected stats
lgen <- list(c(1:15), descrgen$mean, descrgen$sd, descrgen$min,
             descrgen$max, descrgen$skew, descrgen$kurtosis)
lgen <- as.data.frame(lgen)

floor_pct  <- sapply(network_variables, function(x) mean(x == min(x, na.rm = TRUE), na.rm = TRUE) * 100)
ceiling_pct <- sapply(network_variables, function(x) mean(x == max(x, na.rm = TRUE), na.rm = TRUE) * 100)

lgen <- cbind(lgen,floor_pct)
lgen <- cbind(lgen,ceiling_pct)

# Round and rename columns
lgen <- lgen %>% 
  mutate_if(is.numeric, round, digits = 3)

colnames(lgen) <- c("Item", "Mean", "SD", "Min", "Max", "Skew", "Kurtosis", "Floor%", "Ceiling%")

# Display table
kable(lgen, caption = "Table 7. Descriptive statistics for WTM items") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 7. Descriptive statistics for WTM items
	Item	Mean	SD	Min	Max	Skew	Kurtosis	Floor%	Ceiling%
age	1	13.835	1.240	12	18.000	0.303	-0.415	15.744	0.246
AUTO	2	67.583	24.830	0	100.000	-0.574	-0.454	0.861	11.808
HETERO	3	59.513	26.603	0	100.000	-0.148	-0.871	1.353	11.070
TP	4	2.486	0.822	1	5.000	0.298	-0.229	4.182	0.246
SP	5	2.230	0.804	1	4.000	0.237	-0.770	10.578	2.583
RS	6	2.211	0.893	1	4.000	0.290	-0.984	15.867	4.797
CIUS	7	2.427	0.680	1	4.571	0.324	-0.101	0.984	0.123
IDS	8	1.490	0.825	1	5.000	2.105	4.044	52.399	0.492
AB	9	1.484	0.700	1	5.000	2.257	6.099	44.895	0.615
VB	10	1.860	0.899	1	5.000	1.362	1.614	24.477	1.353
VCB	11	1.394	0.705	1	5.000	2.642	8.209	59.410	0.984
ACB	12	1.280	0.595	1	5.000	3.380	13.940	67.651	0.492
AF	13	3.354	0.765	1	4.000	-1.316	1.046	2.091	36.039
AA	14	3.230	0.803	1	4.000	-0.762	-0.554	0.861	34.317
APS	15	3.290	0.745	1	4.000	-0.926	0.056	0.984	33.948

Assessment of multivariate normality assumptions based on Mardia’s tests for skewness and kurtosis, applied to the set of items included in the WTM instrument:

# Run Mardia's tests
skew_result <- mardiaSkew(network_variables)
kurt_result <- mardiaKurtosis(network_variables)

# Round and extract values
mardia_df <- data.frame(
  Test = c("Mardia Skewness", "Mardia Kurtosis"),
  Statistic = round(c(skew_result[2], kurt_result[2]), 3),
  p_value = round(c(skew_result[4], kurt_result[3]), 3)
)
row.names(mardia_df) <- 1:2
# Show table
kable(mardia_df, caption = "Table 8: Mardia's Tests for Multivariate Normality")

Table 8: Mardia’s Tests for Multivariate Normality
Test	Statistic	p_value
Mardia Skewness	5494.117	0
Mardia Kurtosis	48.142	0

2.3. Variables by sex

# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR

item_labels <- c(
  "Age",
  "Autoassertivity",
  "Heteroassertivity",
  "Telepressure",
  "Supervision",
  "Restriction",
  "Problematic use of internet",
  "Suicide ideation",
  "Agression bullying",
  "Victim bullying",
  "Victim ciberbullying",
  "Agressor ciberbullying",
  "Family support",
  "Friends support",
  "Significant Person support"
) 

NV_sex <- NV_sex[NV_sex$sex != "I do not know",]

# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
  item <- items[i]
  label <- item_labels[i]
  
  # Obtener estadísticas por grupo
  stats <- NV_sex %>%
    group_by(sex) %>%
    summarise(
      mean = round(mean(.data[[item]], na.rm = TRUE), 3),
      sd = round(sd(.data[[item]], na.rm = TRUE), 3)
    )
  
  male_val <- paste0(stats$mean[stats$sex == "Male"], " (", stats$sd[stats$sex == "Male"], ")")
  female_val <- paste0(stats$mean[stats$sex == "Female"], " (", stats$sd[stats$sex == "Female"], ")")
  
  # t-test
  t_result <- t.test(as.formula(paste(item, "~ sex")), data = NV_sex)
  t_stat <- round(t_result$statistic, 3)
  df <- round(t_result$parameter, 3)
  p_val <- round(t_result$p.value,3)
    
  # Efecto (Cohen's d)
  d_val <- cohens_d(as.formula(paste(item, "~ sex")), data = NV_sex)$Cohens_d
  d_val <- round(d_val, 3)
  
  # Combinar resultados
  data.frame(
    Item = label,
    Male = male_val,
    Female = female_val,
    T_statistic = paste0(t_stat, " (", df, ")"),
    p_value = p_val,
    Effect_size = d_val
  )
})

# Unir todos los ítems
final_table <- bind_rows(results)

# Mostrar tabla
kable(final_table, caption = "Table 9. Mean (SD) by sex, t-test Results, p-value and Effect Sizes for numeric variables")

Table 9. Mean (SD) by sex, t-test Results, p-value and Effect Sizes for numeric variables
Item	Male	Female	T_statistic	p_value	Effect_size
Age	13.878 (1.251)	13.773 (1.214)	-1.217 (800.604)	0.224	-0.086
Autoassertivity	71.998 (23.962)	63.163 (24.869)	-5.126 (799.408)	0.000	-0.362
Heteroassertivity	65.089 (26.07)	53.8 (25.832)	-6.164 (800.998)	0.000	-0.435
Telepressure	2.41 (0.827)	2.564 (0.807)	2.667 (800.759)	0.008	0.188
Supervision	2.181 (0.799)	2.294 (0.802)	2.008 (800.885)	0.045	0.142
Restriction	2.201 (0.91)	2.232 (0.877)	0.493 (800.316)	0.622	0.035
Problematic use of internet	2.364 (0.638)	2.48 (0.708)	2.437 (791.239)	0.015	0.172
Suicide ideation	1.363 (0.747)	1.597 (0.856)	4.138 (785.192)	0.000	0.292
Agression bullying	1.589 (0.747)	1.355 (0.581)	-4.964 (757.796)	0.000	-0.350
Victim bullying	1.946 (0.959)	1.76 (0.804)	-2.975 (779.264)	0.003	-0.210
Victim ciberbullying	1.403 (0.727)	1.367 (0.644)	-0.729 (791.029)	0.466	-0.051
Agressor ciberbullying	1.319 (0.666)	1.228 (0.459)	-2.264 (714.003)	0.024	-0.160
Family support	3.459 (0.679)	3.268 (0.819)	-3.611 (772.53)	0.000	-0.255
Friends support	3.267 (0.776)	3.21 (0.812)	-1.012 (798.818)	0.312	-0.071
Significant Person support	3.316 (0.727)	3.281 (0.74)	-0.666 (800.519)	0.505	-0.047

2.4. Variables by academic year

# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR

item_labels <- c(
  "Age",
  "Autoassertivity",
  "Heteroassertivity",
  "Telepressure",
  "Supervision",
  "Restriction",
  "Problematic use of internet",
  "Suicide ideation",
  "Agression bullying",
  "Victim bullying",
  "Victim ciberbullying",
  "Agressor ciberbullying",
  "Family support",
  "Friends support",
  "Significant Person support"
) 

# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
  item <- items[i]
  label <- item_labels[i]
  
  # Obtener estadísticas por grupo
  stats <- NV_academic_year %>%
    group_by(academic_year) %>%
    summarise(
      mean = round(mean(.data[[item]], na.rm = TRUE), 3),
      sd = round(sd(.data[[item]], na.rm = TRUE), 3)
    )
  
  male_val <- paste0(stats$mean[stats$academic_year == "First academic period"], " (", stats$sd[stats$academic_year == "First academic period"], ")")
  female_val <- paste0(stats$mean[stats$academic_year == "Second academic period"], " (", stats$sd[stats$academic_year == "Second academic period"], ")")
  
  # t-test
  t_result <- t.test(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)
  t_stat <- round(t_result$statistic, 3)
  df <- round(t_result$parameter, 3)
  p_val <- round(t_result$p.value,3)
    
  # Efecto (Cohen's d)
  d_val <- cohens_d(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)$Cohens_d
  d_val <- round(d_val, 3)
  
  # Combinar resultados
  data.frame(
    Item = label,
    First_Period = male_val,
    Second_Period = female_val,
    T_statistic = paste0(t_stat, " (", df, ")"),
    p_value = p_val,
    Effect_size = d_val
  )
})

# Unir todos los ítems
final_table <- bind_rows(results)

# Mostrar tabla
kable(final_table, caption = "Table 10. Mean (SD) by academic year, t-test Results, p-value and Effect Sizes for numeric variables")

Table 10. Mean (SD) by academic year, t-test Results, p-value and Effect Sizes for numeric variables
Item	First_Period	Second_Period	T_statistic	p_value	Effect_size
Age	13.03 (0.817)	14.911 (0.822)	-32.377 (745.26)	0.000	-2.297
Autoassertivity	68.548 (24.376)	66.293 (25.403)	1.274 (730.568)	0.203	0.091
Heteroassertivity	60.686 (26.028)	57.945 (27.312)	1.444 (727.618)	0.149	0.103
Telepressure	2.434 (0.813)	2.556 (0.83)	-2.101 (739.134)	0.036	-0.149
Supervision	2.439 (0.795)	1.949 (0.729)	9.12 (778.592)	0.000	0.638
Restriction	2.42 (0.874)	1.932 (0.842)	8.037 (761.665)	0.000	0.567
Problematic use of internet	2.392 (0.669)	2.472 (0.692)	-1.66 (733.82)	0.097	-0.118
Suicide ideation	1.526 (0.854)	1.441 (0.782)	1.478 (778.999)	0.140	0.103
Agression bullying	1.453 (0.661)	1.526 (0.748)	-1.435 (693.834)	0.152	-0.104
Victim bullying	1.883 (0.893)	1.83 (0.908)	0.834 (740.582)	0.405	0.059
Victim ciberbullying	1.351 (0.611)	1.452 (0.811)	-1.957 (620.654)	0.051	-0.144
Agressor ciberbullying	1.212 (0.477)	1.371 (0.715)	-3.591 (569.714)	0.000	-0.269
Family support	3.375 (0.761)	3.326 (0.77)	0.907 (742.657)	0.364	0.064
Friends support	3.169 (0.805)	3.31 (0.793)	-2.492 (753.129)	0.013	-0.176
Significant Person support	3.282 (0.731)	3.301 (0.763)	-0.363 (729.875)	0.717	-0.026

3. NETWORK ANALYSIS

3.1. Community Identification

A bootstrap exploratory graph analysis (EGA) was performed to estimate the number, structure, and stability of variable communities (i.e., network dimensions) based on their structural consistency. In this context, a community refers to a cluster of nodes that exhibit strong interconnections through edges.

3.1.1 Community Estimation

The analysis facilitates the identification and delineation of such communities, as well as the allocation of individual items to their corresponding dimensions.

set.seed(123)
communities <- bootEGA(network_variables,
                       model = "glasso",
                       type = "resampling",
                       iter = 500,
                       typicalStructure = TRUE,
                       plot.typicalStructure = TRUE)

Figure 1. Community Estimation

3.1.2. Community Stability

Items were retained within their assigned communities only when their stability coefficients exceeded 0.70. Items falling below this threshold were excluded due to their potential to undermine the structural integrity of the corresponding dimension.

# Structural consistency
communities.dimstab <- dimensionStability(communities)

Figure 2. Community Stability

3.2. Gaussian Graphical Models (GGMs)

3.2.1. Network Estimation

# Definir los grupos (según índices de df_net)
group_list <- list(
  "Parenting Control" = c(1,5,6),
  "Social-Digital Use" = c(2,3,4,7),
  "Bullying" = c(9:12),
  "Social Support" = c(13:15),
  "Suicide"= 8
)

TP_glasso <- qgraph::qgraph(
  cor_auto(network_variables),
  graph        = "glasso",
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 3. Gaussian Graphical Model Estimation

possible_edges <- (ncol(network_variables)*(ncol(network_variables)-1))/2
real_edges <- length(TP_glasso$Edgelist$weight)
density <- round(real_edges/possible_edges*100, digits = 2)

# Print explanation and result
cat("Number of potential edges:", possible_edges, "\n")

## Number of potential edges: 105

cat("Number of real edges:", real_edges, "\n")

## Number of real edges: 19

cat("Network density (percentage of actual edges out of all possible undirected edges):", density, "%\n")

## Network density (percentage of actual edges out of all possible undirected edges): 18.1 %

# Extract edge weights
edge_weights <- TP_glasso$Edgelist$weight

# Calculate summary statistics
mean_weight <- round(mean(edge_weights), digits = 3)
sd_weight <- round(sd(edge_weights), digits = 3)
max_weight <- round(max(edge_weights), digits = 3)
min_weight <- round(min(edge_weights), digits = 3)

# Print explanation and results clearly
cat("Standard Edge Weight \n")

## Standard Edge Weight

cat("The average edge weight was:", mean_weight, "(SD:", sd_weight, ")\n")

## The average edge weight was: 0.178 (SD: 0.298 )

cat("Minimum weight:", min_weight, "| Maximum weight:", max_weight, "\n")

## Minimum weight: -0.245 | Maximum weight: 0.574

# Calculate summary statistics
mean_weight <- round(mean(abs(edge_weights)), digits = 3)
sd_weight <- round(sd(abs(edge_weights)), digits = 3)
max_weight <- round(max(abs(edge_weights)), digits = 3)
min_weight <- round(min(abs(edge_weights)), digits = 3)

# Print explanation and results clearly
cat("\n")

cat("Absolute Values \n")

## Absolute Values

cat("The average edge weight was:", mean_weight, "(SD:", sd_weight, ")\n")

## The average edge weight was: 0.298 (SD: 0.169 )

cat("Minimum weight:", min_weight, "| Maximum weight:", max_weight, "\n")

## Minimum weight: 0.105 | Maximum weight: 0.574

3.2.2. Network characterization

TP_matrix <- as.matrix(network_variables)
invisible(capture.output({SV_mgm <-mgm(
    data = TP_matrix,
    type = rep("g", length(network_variables)),
    level = rep(1, length(network_variables)),
    k = 2,
    verbatim = TRUE,
    warnings = TRUE
  )
}))

#Predice la varianza explicada de cada nodo y otros parametros
SV_predic <- predict(SV_mgm, network_variables, error.continuous = "VarExpl")

TP_glasso_predicted <- qgraph::qgraph(
  cor_auto(network_variables),
  graph        = "glasso",
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 4. Gaussian Graphical Model Characterization

3.3. Centrality Measures

3.3.1. Centality Estimation

centrality_TP <- centralityPlot(TP_glasso, include = 
                                  c("Betweenness","Closeness","Strength","ExpectedInfluence"),
                                orderBy ="Betweenness", scale = "z-scores")

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

Figure 5. Centrality Measures Estimation

3.3.2. Centrality Stability

# 1. Estimar la red con EBICglasso
glasso_net <- estimateNetwork(network_variables, 
                              default = "EBICglasso", 
                              corMethod = "cor_auto",
                              verbose = FALSE)

# 2. Bootstrap de estabilidad (case-dropping)
set.seed(2025)
boot_results <- invisible(bootnet(
  estimateNetwork(network_variables, default = "EBICglasso", corMethod = "cor_auto"),
  type      = "case",
  nBoots    = 200,
  nCores    = 1,        # importante para ver mensajes
  statistics= c("strength", "closeness", "betweenness","expectedInfluence"),  # evita centralidades inestables
  verbose   = FALSE,
  maxErrors = 1000
))

# 3. Plot de estabilidad de centralidad (bootstrap case-dropping)
plot(boot_results,
     statistics = c("strength", "closeness", "betweenness","expectedInfluence"),
     labels = TRUE, order = "sample", legend = TRUE,
     color = "darkblue", line = TRUE,
     main = "Bootstrap Centrality Stability (Case Dropping)")

Figure 6. Centrality Measures Stability

# 4. Coeficientes de estabilidad de centralidad
cs_coeffs <- corStability(boot_results, verbose = FALSE)

# Crear tabla
cs_table <- data.frame(
  Centrality = c("Betweenness", "Closeness", "Strength", "ExpectedInfluence"),
  CS_Coefficient = round(c(cs_coeffs[1], cs_coeffs[2], cs_coeffs[4], cs_coeffs[3]), 3)
)
row.names(cs_table) <- 1:4

# Mostrar como tabla en R Markdown
knitr::kable(cs_table, caption = "Table 10: Centrality Stability Coefficients (CS)",
align = c("r","r","r","r")) %>%
  kableExtra::kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped","hover"))

Table 10: Centrality Stability Coefficients (CS)
Centrality	CS_Coefficient
Betweenness	0.594
Closeness	0.594
Strength	0.750
ExpectedInfluence	0.750

4. NETWORK COMPARISON BY SEX

4.1. Network Comparison by age

TP_male <- NV_sex[NV_sex$sex == "Male",-1]
TP_female <- NV_sex[NV_sex$sex == "Female",-1]

set.seed(1234)
EN_male = estimateNetwork(TP_male, default = "EBICglasso", corMethod = "cor_auto")
EN_female = estimateNetwork(TP_female, default = "EBICglasso", corMethod = "cor_auto")

set.seed(1234)
invisible(capture.output({
  res <- NCT(
    TP_female,
    TP_male,
    p.adjust.methods = "BH",
    binary.data = FALSE,
    test.edges = TRUE,
    edges = "all",
    it = 1000
  )
}))

female_str <- round(res$glstrinv.sep[1], 3)
male_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))

dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))

p_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
    dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
    dplyr::pull(pair),
  collapse = ", "
)

pvalue_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)

In terms of global strength, the female network showed a similary overall connectivity (S = 5.893) compared to the male network (S = 5.745). See Figure 7.

fmt_num <- function(x, d = 3) {
  if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
  if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}

plot(res, what="strength")

Figure 7. P-values for differences in network strength between sexs

However, the network structure invariance test indicated a difference of 0.242 between the female and male networks. This difference was statistically significant (p = 0.049), indicating that the overall network structures differ (Figure 8). Follow-up edgewise tests, with p-values adjusted for multiple comparisons using the Benjamini–Hochberg procedure, showed that the behavior of the edge Item4–Item6 was statistically significant (E = NA, p = NA).

plot(res, what="network")

Figure 8. P-values for differences in network structure between sexs

p_edges

## [1] ""

The Network Comparison Test did not identify any statistically significant differences in edge strength between the UN/NW and OWO networks after applying the Benjamini-Hochberg procedure to correct p-values for multiple comparisons. The significance level was set at the conventional threshold of 0.05.

4.2. Networks by sex

Female’s Network

# --- FEMALE GROUP ---

# Convert data to matrix
TP_female_matrix <- as.matrix(TP_female)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_female <- mgm(
    data = TP_female_matrix,
  type = rep("g", length(TP_female)),
  level = rep(1, length(TP_female)),
  k = 2,
  warnings = TRUE,
  verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_female <- predict(
  TP_mgm_female,
  TP_female,
  error.continuous = "VarExpl"
)

TP_glasso_predicted_female <- qgraph::qgraph(
  cor_auto(TP_female),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 9. Networks of Female Participants

Male’s network

# --- MALE GROUP ---

# Convert data to matrix
TP_male_matrix <- as.matrix(TP_male)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_male <- mgm(
    data = TP_male_matrix,
    type = rep("g", length(TP_male)),
    level = rep(1, length(TP_male)),
    k = 2,
    warnings = TRUE,
    verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_male <- predict(
  TP_mgm_male,
  TP_male,
  error.continuous = "VarExpl"
)

# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_male <- qgraph::qgraph(
  cor_auto(TP_male),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 10. Networks of Male Participants

4.3. Centrality Measures

To evaluate the invariance of the centrality measures network across sexs.

# Compare centralities between networks
network_comparison <- compareCentrality(
  TP_glasso_predicted_female,
  TP_glasso_predicted_male,
  include = "all",
  legendName = "Centrality Measures by sex",
  net1Name = "Female",
  net2Name = "Male"
)
plot(network_comparison)

Figure 11. Centrality comparison between female and male networks

set.seed(1234)
invisible(capture.output({
  nct_test_centrality <- NCT(
    EN_female,
    EN_male,
    it = 1000,
    test.centrality = TRUE,
    p.adjust.methods = "BH",
    centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
  )
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)

# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)

# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15

# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 11. P-values for Centrality Differences (NCT Test)
Node	closeness	betweenness	strength	expectedInfluence
age	0.673	1.000	0.885	0.885
AUTO	0.587	1.000	1.000	1.000
HETERO	0.587	1.000	0.742	1.000
TP	0.587	0.885	1.000	0.587
SP	0.687	1.000	0.587	0.587
RS	0.687	1.000	0.673	1.000
CIUS	0.587	0.587	0.060	0.885
IDS	0.687	1.000	0.673	0.673
AB	0.603	0.901	1.000	0.673
VB	0.741	1.000	0.776	0.587
VCB	0.741	1.000	0.645	0.587
ACB	0.687	1.000	0.587	0.741
AF	0.685	1.000	0.587	0.587
AA	0.776	1.000	0.687	1.000
APS	0.885	1.000	0.587	1.000

5. NETWORK COMPARISON BY ACADEMIC YEAR

5.1. Network Comparison by academic year

TP_first <- NV_academic_year[NV_academic_year$academic_year == "First academic period",-1]
TP_second <- NV_academic_year[NV_academic_year$academic_year == "Second academic period",-1]

set.seed(12)
EN_first = estimateNetwork(TP_first, default = "EBICglasso", corMethod = "cor_auto")
EN_second = estimateNetwork(TP_second, default = "EBICglasso", corMethod = "cor_auto")

set.seed(1234)
invisible(capture.output({
  res <- NCT(
    TP_second,
    TP_first,
    p.adjust.methods = "BH",
    binary.data = FALSE,
    test.edges = TRUE,
    edges = "all",
    it = 1000
  )
}))

second_str <- round(res$glstrinv.sep[1], 3)
first_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))

dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))

p_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
    dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
    dplyr::pull(pair),
  collapse = ", "
)

pvalue_edges <- paste(
  dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)

In terms of global strength, the second network showed a similary overall connectivity (S = 5.243) compared to the first network (S = 5.779). See Figure 7.

fmt_num <- function(x, d = 3) {
  if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
  if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}

plot(res, what="strength")

Figure 7. P-values for differences in network strength between sexs

However, the network structure invariance test indicated a difference of 0.146 between the second and first networks. This difference was statistically significant (p = 0.512), indicating that the overall network structures differ (Figure 8). Follow-up edgewise tests, with p-values adjusted for multiple comparisons using the Benjamini–Hochberg procedure, showed that the behavior of the edge Item4–Item6 was statistically significant (E = NA, p = NA).

plot(res, what="network")

Figure 8. P-values for differences in network structure between sexs

p_edges

## [1] ""

4.2. Networks by sex

second’s Network

# --- second GROUP ---

# Convert data to matrix
TP_second_matrix <- as.matrix(TP_second)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_second <- mgm(
    data = TP_second_matrix,
  type = rep("g", length(TP_second)),
  level = rep(1, length(TP_second)),
  k = 2,
  warnings = TRUE,
  verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_second <- predict(
  TP_mgm_second,
  TP_second,
  error.continuous = "VarExpl"
)

TP_glasso_predicted_second <- qgraph::qgraph(
  cor_auto(TP_second),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 9. Networks of second Participants

first’s network

# --- first GROUP ---

# Convert data to matrix
TP_first_matrix <- as.matrix(TP_first)

# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
  TP_mgm_first <- mgm(
    data = TP_first_matrix,
    type = rep("g", length(TP_first)),
    level = rep(1, length(TP_first)),
    k = 2,
    warnings = TRUE,
    verbatim = TRUE
  )
}))

# Predict explained variance (R²) and other node-level metrics
TP_predic_first <- predict(
  TP_mgm_first,
  TP_first,
  error.continuous = "VarExpl"
)

# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_first <- qgraph::qgraph(
  cor_auto(TP_first),
  graph        = "glasso",
  layout       = TP_glasso_predicted$layout,
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  threshold    = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = SV_predic$errors$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 10. Networks of first Participants

4.3. Centrality Measures

To evaluate the invariance of the centrality measures network across sexs.

# Compare centralities between networks
network_comparison <- compareCentrality(
  TP_glasso_predicted_second,
  TP_glasso_predicted_first,
  include = "all",
  legendName = "Centrality Measures by sex",
  net1Name = "second",
  net2Name = "first"
)
plot(network_comparison)

Figure 11. Centrality comparison between second and first networks

set.seed(1234)
invisible(capture.output({
  nct_test_centrality <- NCT(
    EN_second,
    EN_first,
    it = 1000,
    test.centrality = TRUE,
    p.adjust.methods = "BH",
    centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
  )
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)

# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)

# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15

# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
  kable_styling(full_width = FALSE, position = "center")

Table 11. P-values for Centrality Differences (NCT Test)
Node	closeness	betweenness	strength	expectedInfluence
age	0.455	1.000	0.325	1.000
AUTO	0.325	1.000	0.312	1.000
HETERO	0.325	1.000	1.000	1.000
TP	0.455	0.356	0.280	0.738
SP	0.455	0.710	1.000	1.000
RS	0.455	1.000	1.000	0.312
CIUS	0.280	0.455	0.280	1.000
IDS	0.321	0.781	0.455	0.587
AB	0.587	1.000	1.000	1.000
VB	0.671	1.000	0.455	1.000
VCB	1.000	1.000	1.000	0.582
ACB	1.000	0.240	0.587	0.653
AF	0.280	0.587	0.653	0.280
AA	0.312	1.000	0.672	1.000
APS	0.312	1.000	0.487	0.710

9. BAYESIAN NETWORK

The construction of the Bayesian network (BN) is carried out in two steps. Firstly, the estimation of the directed acyclic graph (DAG) is performed, and secondly, the BN model is fitted and validated with the study dataset.

5.1. DAG Estimation

For the DAG estimation phase, the PC Stable algorithm with no restrictions was employed. To ensure stability in the obtained DAG, a total of 200 bootstrap samples were drawn, and only the edges with a strength greater than 0.85 and a direction greater than 0.5 were retained in the final model.

set.seed(1234)
bootstr <- boot.strength(network_variables, R = 500, algorithm = "pc.stable")

avgnet <- averaged.network(bootstr, threshold = 0.85)

sp <- strength.plot(avgnet, bootstr, shape = "ellipse", render = FALSE)

TP_DAG_qgraph <- qgraph::qgraph(
  sp,
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  nodeNames    = item_labels,
  labels       = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 12. Directed Acyclic Graph (DAG)

5.2. Bayesian Network Validation

To fit and validate the BN model in the second phase of the process, the dataset was subdivided into 10 folds. A routine was implemented in which 90% of the folds were used to train the model, and the remaining 10% were used for testing. This cross-validation routine was repeated until all potential combinations were explore.

#Cross-validation (Check spatial cross-validation: blockCV)
set.seed(1234)
# Split data in 5 sets
kf <- dismo::kfold(nrow(network_variables), k = 10) # k-fold partitioning
models_fit <- c()

for (var in names(network_variables)) {
  kfold_fit <- c()
  aux_base <- data.frame('variable'=var)
  if(is.numeric(network_variables[[var]])){
    for(k in 1:10) {
      test <- network_variables[kf == k, ]
      train <- network_variables[kf != k, ]
      
      training <- bn.fit(avgnet,train)
      pred <- predict(training, node = var, data = test)
      
      z <- data.frame( R2 = R2(pred, test[[var]]),
                       RMSE = RMSE(pred, test[[var]]),
                       MAE = MAE(pred, test[[var]]))
      if(is.na(z$R2)){
        z[is.na(z$R2),] <- 0
      }
      kfold_fit <- rbind(kfold_fit, z) 
    }
    
    z <- apply(kfold_fit, 2, mean)
    z <- cbind(aux_base,as.data.frame(t(z)))
    models_fit <- rbind(models_fit, z)
  }else{
    for(k in 1:10) {
      # Split data in test and train
      test <- network_variables[kf == k, ]
      train <- network_variables[kf != k, ]
      
      training <- bn.fit(avgnet,train)
      pred <- predict(training, node = var, data = test)
      
      validation <- confusionMatrix(as.factor(pred),test[[var]])
      z1 <- as.data.frame(t(validation$overall))
      
      if(length(levels(network_variables[[var]])) > 2){
        z2 <- as.data.frame(t(sapply(as.data.frame(validation$byClass),mean))) 
      }else{
        z2 <- as.data.frame(t(sapply(validation$byClass,mean))) 
      }
      z <- cbind(z1, z2)
      
      kfold_fit <- rbind(kfold_fit, z)
    }
    z <- apply(kfold_fit, 2, mean)
    z <- cbind(aux_base,z)
    models_fit <- rbind(models_fit, z)
  }
}

# Plot the network with R² represented as pie charts on nodes
TP_BN_qgraph <- qgraph::qgraph(
  sp,
  layout       = "spring",
  groups       = group_list,
  #color        = group_colors,
  palette = "colorblind",
  nodeNames    = item_labels,
  labels       = TRUE,
  sampleSize   = nrow(network_variables),
  vsize        = 5,
  esize        = 5,
  asize        = 6,
  border.width = 1,
  border.color = "black",
  unCol        = "black",
  theme        = "colorblind",
  pie = models_fit$R2,
  pieColor = pieColor,
  legend       = TRUE,          # <‑‑ leyenda automática
  legend.cex   = 0.8,           # tamaño texto de la leyenda
  layoutScale  = c(0.8, 0.8),   # encoge ligeramente la red
  rescale      = TRUE
)

Figure 13. Bayesian Network (BN) Model

10. PREDICTAIBILITY CAPACITY

# GGM - General
ggm_general <- SV_predic$errors %>%
  transmute(
    Variable = items,
    R2_gen = round(R2, 3),
    RMSE_gen = round(RMSE, 3)
  )

# GGM - Female
ggm_female <- TP_predic_female$errors %>%
  transmute(
    R2_fem = round(R2, 3),
    RMSE_fem = round(RMSE, 3)
  )

# GGM - Male
ggm_male <- TP_predic_male$errors %>%
  transmute(
    R2_male = round(R2, 3),
    RMSE_male = round(RMSE, 3)
  )

# GGM - First
ggm_first <- TP_predic_first$errors %>%
  transmute(
    R2_first = round(R2, 3),
    RMSE_first = round(RMSE, 3)
  )

# GGM - Second
ggm_second <- TP_predic_second$errors %>%
  transmute(
    R2_second = round(R2, 3),
    RMSE_second = round(RMSE, 3)
  )

# BN model
bn <- models_fit %>%
  transmute(
    R2_bn = ifelse(R2 == 0, "-", round(R2, 3)),
    RMSE_bn = ifelse(RMSE == 0, "-", round(RMSE, 3))
  )

# Combine all into final table
final_table <- bind_cols(
  ggm_general, 
  ggm_female, 
  ggm_male, 
  ggm_first, 
  ggm_second,
  bn
)

# Set column names
names(final_table) <- c(
  "Variable",
  rep(c("R²", "RMSE"), times = 6)  # General, Female, Male, First, Second, BN
)

# Generate final table
final_table %>%
  kable(align = "lcccccccccccc", booktabs = TRUE,
        caption = "Table 12. Explained variance and RMSE of partial-correlation network (GGM) and BN models") %>%
  add_header_above(c(" " = 1, 
                     "General" = 2, 
                     "Female" = 2, 
                     "Male" = 2,
                     "First Period" = 2,
                     "Second Period" = 2,
                     "DAG" = 2)) %>%
  add_header_above(c(" " = 1, 
                     "GGMs" = 10, 
                     "BN" = 2)) %>%
  kable_styling(full_width = FALSE, position = "center")

Table 12. Explained variance and RMSE of partial-correlation network (GGM) and BN models
	GGMs										BN
	General		Female		Male		First Period		Second Period		DAG
Variable	R²	RMSE	R²	RMSE	R²	RMSE	R²	RMSE	R²	RMSE	R²	RMSE
age	0.166	0.912	0.180	0.905	0.146	0.923	0.102	0.946	0.057	0.970	0.116	1.178
AUTO	0.556	0.666	0.533	0.683	0.552	0.668	0.579	0.648	0.536	0.680	0.508	17.332
HETERO	0.483	0.718	0.433	0.752	0.507	0.701	0.473	0.725	0.488	0.715
TP	0.298	0.837	0.382	0.785	0.248	0.866	0.244	0.868	0.385	0.783
SP	0.452	0.740	0.436	0.750	0.461	0.733	0.371	0.792	0.460	0.734	0.415	0.617
RS	0.392	0.779	0.351	0.804	0.431	0.753	0.312	0.829	0.442	0.746
CIUS	0.431	0.754	0.536	0.681	0.335	0.814	0.458	0.736	0.434	0.751	0.414	0.525
IDS	0.373	0.791	0.414	0.765	0.323	0.822	0.412	0.766	0.336	0.814	0.18	0.748
AB	0.555	0.666	0.563	0.660	0.554	0.667	0.508	0.701	0.619	0.617
VB	0.528	0.687	0.601	0.631	0.507	0.701	0.509	0.700	0.561	0.661	0.476	0.652
VCB	0.511	0.699	0.412	0.766	0.578	0.648	0.438	0.749	0.586	0.643
ACB	0.573	0.653	0.451	0.740	0.637	0.602	0.465	0.730	0.666	0.577	0.518	0.406
AF	0.410	0.768	0.353	0.804	0.516	0.695	0.410	0.767	0.454	0.738	0.367	0.606
AA	0.447	0.743	0.512	0.697	0.390	0.780	0.463	0.732	0.408	0.768	0.431	0.607
APS	0.527	0.688	0.545	0.674	0.527	0.687	0.537	0.679	0.544	0.674

Telepressure & Suicide in ESO - Network Analysis

Diego Diaz-Milanes

29/10/2025

0. STUDY CONTEXT AND OBJECTIVES

Reproducibility Statement

Analysis Workflow

1. DATA LOADING & PROCESSING

1.1. Import Packages & Functions

1.2. Data Loading

1.3. Filtering

1.4. Variable Selection

2. DESCRIPTIVE ANALYSIS

2.1. Sociodemographic

2.2. Variables

2.3. Variables by sex

2.4. Variables by academic year

3. NETWORK ANALYSIS

3.1. Community Identification

3.1.1 Community Estimation

3.1.2. Community Stability

3.2. Gaussian Graphical Models (GGMs)

3.2.1. Network Estimation

3.2.2. Network characterization

3.3. Centrality Measures

3.3.1. Centality Estimation

3.3.2. Centrality Stability

4. NETWORK COMPARISON BY SEX

4.1. Network Comparison by age

4.2. Networks by sex

4.3. Centrality Measures

5. NETWORK COMPARISON BY ACADEMIC YEAR

5.1. Network Comparison by academic year

4.2. Networks by sex

4.3. Centrality Measures

9. BAYESIAN NETWORK

5.1. DAG Estimation

5.2. Bayesian Network Validation

10. PREDICTAIBILITY CAPACITY